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ABSTRACT 

We present an algorithm for simulating the equations of ideal magnetohydrodynamics and other 
systems of differential equations on an unstructured set of points represented by sample particles. 
Local, third-order, least-squares, polynomial interpolations (Moving Least Squares interpolations) are 
calculated from the field values of neighboring particles to obtain field values and spatial derivatives 
at the particle position. Field values and particle positions are advanced in time with a second order 
predictor-corrector scheme. The particles move with the fluid, so the time step is not limited by 
the Eulerian Courant-Friedrichs-Lewy condition. Full spatial adaptivity is implemented to ensure the 
particles fill the computational volume, which gives the algorithm substantial flexibility and power. A 
target resolution is specifled for each point in space, with particles being added and deleted as needed to 
meet this target. Particle addition and deletion is based on a local void and clump detection algorithm. 
Dynamic artificial viscosity fields provide stability to the integration. The resulting algorithm provides 
a robust solution for modeling flows that require Lagrangian or adaptive discretizations to resolve. 
This paper derives and documents the Phurbas algorithm as implemented in Phurbas version 1.1. A 
following paper presents the implementation and test problem results. 
Subject headings: Magnetohydrodynamics (MHD), Methods: numerical, Hydrodynamics 



1. INTRODUCTION 
1.1. Context 

Our understanding of many astrophysical systems re- 
lies on the simulation of magnetized plasmas. As a re- 
sult, much effort has been made to develop tools to effi- 
ciently perform high-fidelity simulations of them. Some 
of these tools have found broad application in other fields 
of physics and engineering as well. 

Many early methods for solving the equations of mag- 
netohydrodynamics (MHD) were based on fixed grids. 
Discretizing the equations of hydrodynamics or MHD on 
a fixed grid leads to an Eulerian method, or a method 
written in terms of Eulerian derivatives. Popular pub- 
licly available codes with methods based on point values 
such as the Pencil Code F] ( (Brandenburg fc Dobler.2002| ) 
and finite volumes, sucn as ZEUS dHayes et al.||'iOO(j|)' 
FLASH ( [Fryxell et al]|2000[ ), or Athena ( |stone et al. 
20081 use such methods! Eulerian methods share the 



common property that the discretized form of the gov- 
erning equations is not Galilean invariant. Though they 
still converge to the correct solution, this does lead to 
two limitations at any finite resolution. First, the ex- 
plicit integration time step constraint from the Courant- 
Friedrichs-Lewy (CFL) condition depends on both the 
signal speed and the fiow velocity relative to the grid, 
not just the signal speed. Second, the numerical diffu- 
sion of the scheme, usually highly nonlinear, also depends 
on the fiow velocity relative to the grid. 
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A fixed grid approach thus has disadvantages partic- 
ularly where there are high-velocity bulk flows, collaps- 
ing flows, or flows that generate localized fine str ucture. 
For the latter c ases, adaptive mesh refinement (Berger 
& Oliger 1984) has been a successful approach, 
method 



uses refined meshes to al- 
How- 



while still Eulerian, 
low the spatial and temporal resolution to vary, 
ever, for problems with significant bulk flows, it is of 
no help, as the same problems of time step limitation 
and numerical diffusion apply as with uniform grids. A 
numerical viscosity dependent on the bulk flow can be 
significant, because the growth of instabilities from a 
marginally resolved mode in a method lacking Galilean 
invariance will depend on the bulk velocity of the fiow 
across this grid. The effects of this can be seen, for 



example, in Chiang (20081 and Johansen et al. (2009). 
To circumvent the time step liniit m disks treated with 
cylindrical or spherical coordinates or in a shearing-sheet 
approximation where the bulk flow is largely Keplerian 
and aligned with the grid, it is possible to ad d a separate 
transport step to the method ( Masset|[2000 ) . While this 
extra transport step improves the problems wi th numer- 
ical diffusion, it does not fully cure the issue (IJohansen 



et al. 112009 Stone & Gardiner 



2010) 



^ See http://wMW.nordita.org/software/pencil-code/ 



lb escape these limits, it is necessary to move to a 
method formulated in terms of Lagrangian (also known 
as covariant, comoving, convective, advective, substan- 
tive, or material) derivatives]^ In contrast to Eule- 
rian formulations, Lagrangian methods have three ad- 
vantages. Foremost, for problems with significant bulk 
flows, a purely Lagrangian formulation has a significantly 
less stringent time step constraint from the signal speed 

* Methods that solve Eulerian problems in a local frame chosen 
to be comoving with the fluid in a locally average sense also share 
in some of the advantages of this formulation. 
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(the CFL condition). This is because the time step in an 
Eulerian method depends on the maximum of the signal 
speed and the flow speed, whereas in a purely Lagrangian 
method the time step depends only on the local signal 
speed. Relaxing this constraint becomes particularly im- 
portant in the case of an extended disk with supersonic 
differential rotation, where in an Eulerian formulation 
the quickly orbiting inner regions constrain the time step 
severely. A second advantage of a Lagrangian method 
lies in the Galilean invariance of the inevitable effects 
of numerical diffusion. Though Galilean invaria nce itself 



can formally be achieved in an Eulerian method ( Springel 
2010 Robertson et al. 2010), a Lagrangian formulation 



can reduce the dittusivity further because it uses fewer 
time steps. Finally, Lagrangian methods naturally focus 
resolution into regions of fluid concentration, which are 
often, though not always, the regions of greatest interest. 
(We note that the adaptive, Lagrangian method we de- 
scribe here can also focus resolution to other, arbitrary 
regions of interest.) 

It is possible to write a comoving discretization in two 
ways. First, one can discretize the governing equations 
directly in terms of Lagrangian time derivatives. Second, 
one can discretize in terms of partial time derivatives 
around moving interfaces. Historically, the most popular 
approach has been the first, particularly when used to 
build a meshless method. Recently, the second has been 
used, with techniques based on a moving imstructured 
mesh with mesh reconnection. 

One of the earliest and most popular meshl ess schemes 



is Smoothed Particle Hydro dynamics (SPH; Lucy 



1977 



Gingold & Monaghan 19771. SPH quickly gained pop 
ularity as the advantages m numerical diffusion, local 
resolutio n scales, and local time st ep advantages were 
realized ( Steinmetz fc Mueller|199^ ). However, the basic 
SPH algorithm has many short cornings. The foremost 
and most fundamental is the lack of discrete, zeroth- 
order consistency in the SPH representation of a func- 
tion. SPH interpolation fails to reproduce even a con- 
stant function. The importance of this consistency prop- 
ert y in general mesh less schemes has been pointed out 
by Liu et al. (1995 ). This insi g ht has been app li ed to 



analysis of SPH"^ ypiIts| (|1999[); [Liu et al. 
fc Matthiesl ( |2004[ ) aH ^umlEi et al.| ( |200(j' 



(2003[); [Fries 
among oth- 
ers. I'hey tind that the lack of zeroth-order consistency 
can cause substantial gradient and value errors that do 
not converge with increased particle number alone. The 
inability of SPH to effectively model subsonic t urbulence 
has been blame d on this lack of consistency by [Bauer fc| 
Springel (2011), though the behavior in this regime de- 
pends strongly on, and can be significantly improved by 
using a mo re rnodern formulation of the SPH artificial 
viscosity (Price[ [20lT ). 

Resolution in sPH is further limited by constant par- 
ticle masses. Soi ne attempts at adaptive parti cle masses 
have been made ( Kitsionas fc Whitworth|2002 ) but these 
suffer from difficulties in specifying a well-posed scheme. 
SPH in general handles differing particle masses poorly, 
as the pairwise interparticle interactions allow heavy par- 
ticles to penetrate though the fluid in a nonphysical man- 
Similarly, the spatial resolution in SPH is locally 



ner 



isotropic, even when the particle and mass distribution 
is anisotropic. Attempting to r elax this constraint leads 
to the adaptive SPH scheme of Shapiro et al. ( 1996 ) and 



Owen et al.[ ( [1998[ ). 

A grid that is both Lagrangian and has logically Carte- 
sian structure is a simple choice, and a logically Carte- 
sian moving (Lagrangian) mesh has also been used to 



attempt t o minimize num erical diffusion (Norman et al 
198(3. Fiedler fc Mouschovias„1992) . .Gncdm, (1995] and' 
ren 1998) used a moving, logically Cartesian mesh to 
provide adaptivity in collapsing flows. However, this ap- 
proach falls victim to several limits. In many flows the 
cells eventually become long and thin, leading to large 
errors. Also, the grid cannot follow rotation or turbulent 
flows as it becomes tangled. 

Unstructured, moving mesh methods with mesh recon- 
nection have re cently bee n intro duced in astrophysics . 
The methods o f 'Spring ej ([2 010|, [Pakmor et alT (1 2011 ) 
Duffell fc MacFa dycn (2011 ' 
are finite volume methodsl 



and Gaburov et al 
^sed on Voronoi tessellations 
The mesh is defined by the Voronoi tessellation of a set 
of points that move approximately with the mean mo- 
tion of the fluid in the cell (though formally any mo- 
tion can be chosen). These methods can be described as 
Lagrangian though they calculate inter-cell fluxes with 
Eulerian Riemann problems stated in a locally comov- 
ing frame. The connectivity of the mesh is dictated by 
the Voronoi neighbor relation. Fluxes between cells are 
calc ulated across the mov ing cell faces 
( [201l|_ 
I'he method of lu 



and Pakmor et al 
method. 



Springel[ ( [2010 1 
describe a Galilean invariant 



uffeU & MacFadycn 
not fully Galilean invariant, but this is due to t 



(|2011| is 

le formu- 



lation chosen for the slightly more complicated relativis- 
tic hydrodynamic equations. Both methods use an ap- 
proximately comoving formulation in a significant sense. 

This paper describes an adaptive, Lagrangian, mesh- 
less, collocation scheme for MHD or similar sets of equa- 
tions based on a point (not finite volume or mass) dis- 
cretization. In what follows, we refer to the discretization 
points as particles, following the historical usage. How- 
ever, these discretization points do not in any sense rep- 
resent identifiable masses or volumes of the fluid. They 
are simply moving points sampling continuous fleld vari- 
ables. 

In the next subsection we discuss prior work on re- 
lated methods to solve the MHD equations. We then 
describe our algorithm, starting with an overview (§ [2| 
and then discussing speciflc numerical aspects, such as 
the modeling of the function and the time update (§ [3|), 
adaptive addition and deletion of particles (§[4|), explicit 
time step limits (§[5|), and magnetic divergence correction 



(§ 6.2). Finally we draw these together with a summary 



of the algorithm (jj| 7|. In the next paper of this series 
( McNally et al.[[20l2} hereafter Paper II) we present im- 
plementation details and present the results of a suite of 
gas dynamical and MHD tests of the algorithm. 

1.2. Prior MHD Methods 

Several attempts have been made to design an SPH- 
type sch eme for MHD. The most suc c essful and recent 
work by [Price fc Monaghan[ ([2004a|b[ [2005[ ), and [Price 



( 2010| resulted in an SPH MHD based on a form of the 
IvlHD equations that is consistent wi th V ■ B and 
a set of artiflcial dissipation terms. [Rosswog fc Price] 
(20071 developed a variation based on representing the 



magnetic fields though Euler angles, which allows a guar- 
anteed V • J5 = at the cost of disallowing tangled field 
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geometries (|Brandenburg|2010 ), seve rely limiting its ap- 
plicability. Dolag &: Stasyszyn ( 2009 1 implement an SPH 
MHD in GADGE'i'-a, without any constraint on V • B, 
but subtracting the numerical contribution of V • B to 
the m omentum equation. We refer the reader to |Price| 
(2012) for a further overview of the attempts to design 
aFKPH MHD method. 

Unfortunately all these SPH MHD methods suffer from 
the fundamental drawback of SPH, that the SPH inter- 
polant does not have a zeroth-order consistency prop- 
erty. This zeroth-order inconsistency means that for a 
disordered set of SPH particles, a constant function can- 
not be reproduced by the SPH representation of that 
function. As the SPH representation of even a constant 
function has significant positive and negative errors, it 
also has significantly non-zero derivatives. These errors 
make formulating an SPH MHD difficult. Modifications 
of SPH to solve or work around the zeroth-order consis- 



tenc y problem have bee n proposed. B0rve et al. (20011 
and B0rve et al. (2006) developed an extension to SPH 



using a remappmg strategy to increase the accuracy of 
SPH estimates through regularizing the particle distribu- 
tion, and appl ied it to MHD sho cks. For hydrodynamics, 



Morris] ( |1996[) and Abel ( 2011[ ) have proposed working 
around the effects of the zeroth-order consistency prob- 
lem for pressure forces only, with an alternative deriva- 
tion of the SPH pressure force. This comes at the price of 
sacrificing the local momentum conservation enjoyed by 
the classical formulation. This also only treats the prob- 
lem of spurious pressure forces arising from the zeroth- 
order inconsistency, and does not lead to a consistent 
interpolation of the pressure field or other fields. 

It is also possible to construct a SPH MHD scheme us- 
ing a Godunov approach. Godunov SPH was originally 
proposed for hydrodynamics, using Riemann problems 
to solve for the particle interactions. Godunov SPH uses 
SPH interpola. tion for densit y (see Eg. 6 and Eq. 21 of 



Inutsuka 2002 and Eq. 29 of Iwasaki & Inutsuka 2011) 



^ A Godunov SPH MHD implementation using Powell- 
type sourc e terms and a tensile corr ection was imple- 
mented by Iwasaki & Inutsuka (2011). They point out 
that all SPH- based MHU schernes that avoid tensile in- 
stability do not exactly conserve momentum , energy, or 
both. Similarly, Gaburov & Nitadori (2011) constructed 
an SPH-like scheme (a weighted particle method) with a 
consistent second order accurate formulation for deriva- 
tives, coupling this with a pairwise Riemann-solver based 
interaction between particles to y ield an MHD s c heme . 
A Galilean invariant form of the Dedner et al. (2002) 
hyperbolic-parabolic cleaning scheme was used to han- 
dle V-B errors. 

However, these SPH-based methods again suffer from 
the zeroth order inconsistency of the SPH interpolant, 
even though methods with a renormalized first derivative 
estimate have a consistent first derivative. This means 
that SPH interpolated fields (such as the density values) 
have significant noise. To reduce the amplitude of the 
noise it is necessary to increase the number of neighbors 
used in the kernel, which greatly increases the computa- 
tional cost. This means that rigorous convergence stud- 
ies, even in smooth flow, are not feasible with methods 

An earlier usage of Ri emann solvers coupled with SPH is given 



An earlier usage ot Kis 
by |Parshikov et al.|f2000^ 



based on SPH-type estimates. In addition, SPH Rie- 
mann methods suffer a higher computational cost in com- 
parison to moving unstructured mesh Godunov schemes, 
because of the requirement of a much higher number of 
Riemann p roblem solutio n s per particle. 



iDuffell fc Ma cFadyen 
their 



(2011 



implemented a MHD 
Voronol tessellation method, using a 

but 



scheme in 

Dedner type hyperbolic divergence cleaning method 
found it difficult to manage \/-B errors w hen the mesh 



topology changes. Pakmor et al. (2011) used a very 



similar approach, with ap parently much g r eater success 
in managing W B errors. Gaburov et al. (2012) add a 
source term to the induction equation to restore Galilean 
invariance if V-J5 7^ and claim this greatly improves 
stability. 

The method we descri be here was inspired by the Gra- 
dient Particle Method of |Maron fc Howes] (2003), but re- 
moves the underlying instability present m that method 
(descri bed in Appendix K). A m ethod particularly sim- 
ilar to iMaron fc Howes ~! 2003 ) , but limited to hydro- 



dy namic s using a mo vmg-least-squares fit was proposed 



by Dilts (1999 2000). Numerous related methods have 



been described m the literature on meshfree or meshless 
methods. The most closely related met hod is the Finite 
Pointset Method (FPM) described by ]Kuhnert] (]1999 



2002 ) , which is not to be confused with either the simi 



larly named Finite Point Method of|Onate et al. 1996 



or the equally similarly named Finite Particle Method o: 



Liu et al.| |200fp] FPM has limited adaptivity, is first or- 
der, and uses an upwinded formulation for hydrodynam- 
ics. Similar to the method we describe, it is meshless, 
Lagrangian, has particle addition and deletion, and uses 
moving least squares interpolation. 

2. ALGORITHM 

For specificity, we focus on using our method to solve 
the equations of MHD. These can be expressed using 
Lagrangian time derivatives Dt, as 



P ^£jab£acd{dcBd)Bb + G 



J' 



(1) 

DtB, = B,d,V, - B,d,V^, (2) 



Dt<J = -{a + P)d^V^ 

Dtp = -p^^V^, 



(3) 
(4) 



where V is the velocity, B is the magnetic field, a is 
the internal energy volume density, P is the pressure, 
p is the density, Gj is a vector component of a body 
force, and the Einstein summation rule is assumed. We 
note that Phurbas is relatively insensitive to the exact 
form of the equations solved and variables chosen. For 
example, energy variables other than the internal energy 
per volume could be used. In Appendix ]C] we give the 
second time derivatives of these equations tor use in the 
time update. These equations require the addition of an 
equation of state, such as a gamma-law P = (7 — l)cr, 
though the equation of state is arbitrary. 

The MHD equations (ll])-(]4]) are solved on an adaptive 
set of particles, each particle carrying values for the field 

^ The authors are of the opinion that enough numerical schemes 
have been named FPM, and as the names are getting confusing the 
practice should cease. 
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variables p, V, B, and a. Particles move in the frame of 
the fluid with the local fluid velocity V. Field variables 
evolve in the frame of the particle, so the evolution equa- 
tions are most naturally expressed using the Lagrangian 
form for the time derivatives in the MHD equations. 

The equations of MHD as stated are ill suited to the nu- 
merical scheme we will use. For the discretization used in 
Phurbas, we require a system of equations in which short 
wavelength perturbations decay. Appendix [A] demon- 
strates this for a model advection-diffusion equation. To 
ensure decay of such perturbations, we introduce artifi- 
cial dissipation terms to the analytic form of the equa- 
tions before discretizing. These modifications are in the 
form of a bulk viscosity, and mass and thermal diffusions. 
Formally, it is this modified version of the MHD equa- 
tions from which Phurbas computes approximate numer- 
ical solutions. The MHD equations, reiterated with the 
addition of the stabilizing terms, and associated fields 
are: 



DtB, 
Dtu 



^djP + p ^£jab£acd{dcBd)Bb 

+ G,+^,{{C, + Cl)^^V,), 

+ H^pdiiCsdA), 
P 



Dtp=^pd,V, + Hpd,{Csd,p), 

DtCl = di{KQdiCl) + — ACmax " 
Tl 



1 



DtCs = d^{K^dtCs) H Ss —Cs, 

Ts+ Tg 



(5) 
(6) 

(7) 
(8) 
(9) 

(10) 



where A is the Nyquist length, and Cmax = \/cfT^ is 
the maximum signal speed, where Cg is the sound speed 
and Va is the Alfven speed. Bulk viscosity fields Q and 
Cs are introduced to handle the general flow (Q), and 
shocks and other discontinuities (Cs)- The action and 
parameterization of these fields are described in § |6.1 
H„ and Hp are constants that specify the strength o: 
mass and thermal conductivities in continuity and energy 
equations, while is the term representing diffusion of 



6.2 



magnetic divergence defined in 

To evolve the field variables in time, we evaluate Equa- 
tions (|l}-(|4]) for the time derivatives, requiring values 
for the neldvariables and their spatial derivatives at the 
position of each target particle. We obtain this infor- 
mation by fitting a third-order, three-dimensional (3D) 
polynomial to the set of values carried by the neighbor- 
ing particles, using the procedure described in § [Sj The 
resulting polynomial coefficients allow us to compute the 
field value and its first, second, and third derivatives at 
the position of the particle, enabling evaluation of the 
Lagrangian time derivatives. Those in turn are used to 
update the field variables with a predictor-corrector time 
step scheme described in 



3.3 



A particle-based algorithm such as this one has a dy- 
namically evolving spatial resolution. It turns out to be 
central to the stability and accuracy of the method that 
the particle distribution not have voids within which the 
fields cannot be accurately fit. We create and delete par- 
ticles as necessary to eliminate such voids, while avoiding 



particle clumps. This further allows us to adaptively sat- 
isfy any user-specified physical resolution requirement, as 
well as to eliminate unnecessary particles (§[4]). We force 
the resolution to always exceed a spatially and tempo- 
rally variable target resolution X{x,t). Effectively, the 
particles can represent the field variables in the same 
manner as a grid with effective resolution A at each point. 
The resolution requirement can be specified depending 
on the physics requirements of the problem at hand, so 
long as it remains reasonably smooth. 

The Phurbas discretization is based on point val- 
ues, not finite volumes or finite masses. As such, 
the discretization used to calculate spatial derivatives 
and advance the solution in time does not define a 
value for volume-integrated quantities, including volume- 
integrated, conserved quantities. To define these quan- 
tities, another discretization would need to be added to 
obtain a multidimensional quadrature from the unstruc- 
tured set of samples. For example, Voronoi cell volumes 
could be used to calculate a nearest-grid-point interpola- 
tion for a Riemann-sum approximation to a volume inte- 
gral. Alternatively, using a point density approximated 
from the number of neighboring points in some small ra- 
dius as a weighting, a Monte-Carlo type volume integral 
approximation could be used. 

As the magnetic field evolves, discretization error gen- 
erates spurious magnetic divergence V B. By dropping 
the physically vanishing term —V B when deriving the 
Lagrangian induction equation from the usual Eulerian 
form expressed with a partial time derivative, we have 
made any V-B present in the field into a passively ad- 
vected scalar. A consequence of this choice of the canoni- 
cal Lagrangian form is that our MHD equations, by omit- 
ting a term that is physically zero, are precisely the same 
as a form that is claimed in other works to include an ex- 
tra source term: the same result h as been proposed with 
a source term by Janhunen (2000), and derived from the 
relativistic form of energy-momentum conservation and 
relativistic electromagnetic theory by Dellar (2001). In 
the latter paper, it is shown to be the Galilean invariant 
momentum and energy conserving form for the MHD 
equations in the case when V-B is present. We note 
that V-B of nonphysical origin is easier to numerically 
handle in the Lagrangian form of the MHD equations, 
as the presence of V-B errors does not feed back into 
violations of energy and momentum conservation by it- 
self. Unhkc in SPH MHD iPricCj ( ,2012) , we do not require 
source terms in the momentum and energy equations, as 
our discretization does not suffer from the tensile insta- 
bility as in SPH. The diffusive correction described in 



6.2 adds a term to the magnetic equation that diffuses 



3. TIME EVOLUTION 



We evolve the field variables forward in time by evalu- 
ating the MHD equations ([l])-(|4]) at the position of each 
particle. We do this by constructing a local approxima- 
tion of the field variables at that position, derived from 
a spatial fit to the valu es of the field variables on neigh- 
boring particles (§ 3.1). This allows us to compute the 
values and spatial derivatives of the field variables at 
the position of the particle. We choose for the form of 
the continuous approximation a 3D, third-order, polyno- 
mial. We further develop a system of particle weights 
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that enhances the accuracy of the fit (§ 3.2). Once we 



have evaluated the Lagrangian time derivatives from the 
MHD equations, the field variables and particle positions 
ar e up dated in time with a predictor-corrector method 



(§ 3.3) 



3.1. Moving Least Squares Procedure 

Phurbas uses a moving least squares fitting procedure 
in two versions. First, to approximate derivatives of the 
dynamical field variables on the right hand sides of the 
governing equations, moving least squares interpolants 
are used. These are polynomial approximations that are 
forced though (interpolate) the central particle value. 
Second, Phurbas uses moving least squares fits to ini- 
tialize newly created particles. These are polynomial ap- 
proximations that are not forced though any particle val- 
ues, and hence provide a smooth approximation to the 
field values where no particle currently exists . 



In the language of Lancaster & Salkauskas ( 1981 ) the 
first version is an Interpolating Moving Least Squares 
procedure, and the secoi id is a Moving Least Squa r es pro - 
cedure. In addition to Lancaster & Salkauskas (1981[, 
discussion at length c an be found , fo r example, in [Be 



lytsch ko et al. (1996), Dilts (1999), or Fries & Matt 



hies 



I'tar t he purposes o f Phur bas, we can follow the descrip- 
tion by Liszka et al. (19961, which leads to what has be- 
come known as a Gene ralized Finite Difference Method 
(iLiszka & Orkisz||1980'). P hurbas uses both the Nodal 
Approximation described in Liszka et al. ( 1996 section 
2.1.1) and the Pointwise Approximation m Liszka et al. 
(1996 section 2.1.2). We briefly expand on those de- 
scriptions here for clarity. 

In one dimension, we start with a function that we wish 
to discretize, defined at some set of points g — g{x.). 
To approximate the function at a point Xq, we select 
a set of nearby points {x^}, and then write the series 
approximation using n polynomial terms Pi(x), 



9(x) = ^aiPi(x-Xo) 



(11) 



so that g(x) « .9(x). If the functions Pi are selected as 
polynomials. 



p(x) = [l,x,y,z,x'^,y'^ 



2 3 



(12) 



then this approximation is a Taylor series about xq, as 
in Equation (11) we have shifted the polynomials p to 
be centered on xq. The coefficients Oi are then the value 
and derivatives (multiplied by a Taylor series coefficient) 
of the approximation q{x). 

If the function g is defined at Xg we can reduce the 
approximation to a special case, c alled the Nodal Ap- 



proximation by Liszka et al. (1996 section 2.1.1). If we 



we seek a solution in the least-squares sense. That is, the 
solution for should minimize the quadratic form 



where is a weight function described below. We can 
rewrite this set of equations in matrix form by defining 

g"^ = [5(xi) -5(xo), 5(^2) -5(xo), 5(^3) -5(xo),---] , 



(15) 



Pl(Xi-Xo) P2(Xl-Xo) ... p„(Xi-Xo) 
Pl(x2-Xo) ;72(x2-Xo) ... Pm(x2-Xo) 



Pl(x„-Xo) P2(x„-Xo) ... p,„(x„-Xo) 



(16) 



and 



W 



VK(xi - xo) 

I^(X2 - Xo) 














M^(x„ - Xo) 



Then Equation ( 14 ) can be written 

J = (Pa-g)'^W(Pa-g). 



If we define 



A ^ P^WP 

B = P^W 



(17) 

(18) 

(19) 
(20) 



then minimizing J as in Belytschko et al. ( 1996 ) we ob- 
tain 



and 



a= A-^Bg 



dq dq dq 2 
(9a; ' dy' dz' dx"^ 



(21) 



(22) 



gives the derivatives of the interpolating moving least 
squares approximation. 

The second form, the Pointwise Approximation, is de- 
fined everywhere, not just where there is a particle. This 
form is used when adding new particles. The coefficient 
ao is left free, the approximation is Equation (11), and 



fix the coefficient ao = 3(xoj = g(xo), then only the co- j^^^ ^t an arbitrary point x the approximation ^I^l'ds 



efficients ai for z > 1 need to be determined, and the 
approximation becomes 



g(x) = 5r(xo) + ^ ajpi(x - Xq) 



(13) 



We choose to use a number of neighboring points n 
greater than the number of und etermined polynomial co- 
efficients TO. As Equation ( 13 1 is then overdetermined. 



= [.9(xo),5(xi),g(x2),g(x3), . 

Po(xo - x) pi(xo - x) 
Po(xi - x) pi(xi - x) 

Po(x„ -x) pi(x„ -x) 



(23) 



p™(xo - x) 
p™(xi - x) 



(24) 
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and 



W = 



VK(xo - x) 

W{yii ~ x) 



The vector of coefficients then yields 



^ ^ dq dq dq ^^^1 



(25) 



(26) 



The coefficients vector has either 19 or 20 coefficients, 
depending on the approximation, so the solution requires 
inversion of either a 19 x 19 or 20 x 20 matrix A. We use 
an L U decomposition and back substitution proce dur e 



(e.g. [Press et al.]|1992[ p. 32) to solve Equation 
The derived polynomial coefficients yield the values of 
the field variables and their derivatives of first, second 
and third order, from which we construct the first- and 
second-order time derivatives of the field variables. 

We have experimentally found that the number of par- 
ticles included in the evaluation sums for the matrix coef- 
ficients should comfortably exceed the number of terms in 
the polynomial. The choice of how many particles to in- 
clude is based on a compromise between lack of statistical 
significance and computational impracticality. The ra- 
dius Tf oi the sphere encompassing the particles included 
should also be large enough to justify a third-degree in- 
terpolation, about twice the characteristic inter-particle 
separation A. For a uniform particle density of one par- 
ticle within each volume A^, a sphere with rf — 2X en- 
closes ~ 34 particles. However, this leaves little room to 
account for non-uniform particle densities. A sphere with 
rf — 3X encloses ~ 110 particles, which weighs on the 
cost of calculating the matrix coefficients. A radius as 
large as this also invites higher-order structure to erode 
the interpolation. In the end we choose to use 



2.3A, 



(27) 



corresponding to ^ 51 particles. We have not yet de- 
rived a rigorous lower bound to the required number of 
neighbors to use. For computational cost, we find that 
a third-order fit has computational cost comparable to 
the other operations that occur in the time step, while a 
fourth-order fit is substantially more expensive. We term 
the sphere of radius r f around a fit center the neighbor 
sphere. 

The target resolution X{x,t) is a property of the loca- 
tion of the fit center, which may or may not be centered 
at the location of a particle. In practice, as the interpo- 
lation is centered at the location of an existing particle, 
then A for the interpolation is taken as the A of that par- 
ticle. If a fit is being used to generate new values for the 
addition of a particle, the A of the particle that triggered 
the creation of the new particle is used. 

3.2. Weights 

In the moving least-squares procedure, each neighbor 
particle j has a weight Wj. There is significant free- 
dom to choose the form of the weights Wj , and no rigor- 
ous theoretical framework exists under which an optimal 
choice can be made. As a first approximation, we choose 



weights that emphasize particles close to the center of 
the neighbor sphere, so that the local least squares ap- 
proximation varies more smoothly as the target position 
is changed. The weighting function is a piecewise linear 
function of the distance of particle j from the fit center 



if < rj < r„ 



(P^) ifr. <r, <r^. 



(28) 



where r^, = |A. Wj has a nonzero value at the edge of 
the neighbor sphere so as to not exclude any particles 
from the fit. 

3.3. Time Update 

The field variables are evolved in time with a Her- 
mite predictor-corrector scheme based on the first- and 
second-order Lagrangian time derivatives. A derivation 
of the scheme is presented in Appendix |B] as it has not 
previously been described in the literature. The interpo- 
lation procedure in § 3.1 for time step i -I- 1 is done on 
the predicted values qp^i computed in time step i, yielding 
the time derivative values Dtqp^i and Duqp,i needed for 
the correction in time step i, as well as for the prediction 
of time step i + 1. 

We begin by extrapolating forward from time ti to time 
ii+i, over the time interval At = t^+i — ti, to make a 
prediction 



Dtqp,^At + lDuqp,^At^ , (29) 



based on a Taylor series expansion around qp^i, using 
the corrected value from the previous time step Qci- We 
then evaluate the time derivatives by interpolation on 
the predicted fields g^.i+i at time t^+i, and correct the 
prediction to derive the corrected value at ti+i. 



qc,i+i ^qcA + ^{Dtqp,i + Dtqp^i+i)At 
[DttqpA - Dttqp,i+i)At'^. 



(30) 
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The particle positions x are evolved using third-order 
time information DutX = Dtt V as well. This allows us 
to use a third-order predictor of the form 



X. 



'p,i+l 



^Xc,i + Vc.iAt 

^DtVp^.At^ + Id 



(31) 



2 6 
and to correct it to the final value 
1 



uVp,,At^ 



+ M.^ + Vc,^+l)At 



1 



10 
1 

120 



{DtVp,, - DtVp,,+,)At^ 



(32) 



{DttVp,^ + DuVp^^+i)At'' 



4. REGULARIZING THE PARTICLE DISTRIBUTION 

Our algorithm relies on discretization over Lagrangian 
sample points. These points are not arrayed on a grid, 
nor are they connected by mesh edges as in the AREPO 
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code ( Springel|2010 ), so this is a meshless method. Dur- 
ing evolution, we require that the particles should main- 
tain a distribution such that there are no voids larger 
than the target local resolution X{x,t) and no excessive 
point concentrations within the scale A. The requirement 
of no voids is introduced to ensure that every fit sphere 
of radius rj has enough points to perform the moving 
least-squares procedure and that the fit spheres overlap 
sufficiently so that the set of fit spheres covers the entire 
simulation volume. The requirement of no point concen- 
trations requires the removal of excess points. We im- 
plement these requirements by adding and deleting par- 
ticles as needed. Satisfying these requirements confers 
the great benefit of making the code fully adaptive, since 
the user can dynamically choose the function X{x,t) as 
required by the physics of the problem, so long as it is 
reasonably smooth in space and time. 

The addition and deletion algorithm begins with the 
assembly of all neighbors i within the neighbor sphere of 
radius rf around a particle j, along with their associated 
target resolutions Ai. Voids within the neighb or s phere 
are identified using the method described in § |4.1| Any 
voids identified are reported as candidates for particle 
creation. Conversely, if a partic le 7 has a mutual nearest 



neighbor that is too close (see § |4.2[ ) , one of the two par- 
ticles is deleted. Duplicate voids and clumps are pruned 
from the global list prior to the particle creation and 
deletion described in § 14.31 



4.1. Voids 

To check for a void at a point in space with position 
X, we identify the nearest particle i, which is located at 
position Xi and has a resolution scale A^. The distance 
between x and the particle normalized by the resolution 
scale is then 



I X Xi 

A,, 



(33) 



As a resolution condition, we then choose the condition 
that if 



•^void ^ ^voidi 



(34) 



for a constant Cyoid, the space around x is indeed too 
sparsely populated, indicating a need for particle addi- 
tion. 

To heuristically derive Cyoid, we consider an arrange- 
ment of particles on the hexagonal lattice representing 
the tightest possible packing of spheres centered on the 
particles. If the particle density is one particle per vol- 
ume A'^, then the spheres' centers will be separated by 
a distance dp — 2^/^. This represents the most efficient 
possible filling of the region with particles. Since any 
real, fluctuating, particle distribution will require more 
particles to fully resolve the field, we set 



^void — 0.73fir) 



(35) 



so that with a disordered particle set we sample more 
density than would be required with the ideal ordered 
particle set. 

To identify unique voids, we first identify the most 
egregious void within the neighbor sphere of each par- 
ticle, and then check to see if that void violates the con- 



as described below in Section [4.3[ We begin by examin- 
ing the space in the vicinity of existing particles. We 
construct a 3D cubic grid with side length 2r/ contain- 
ing 9x9x9 grid points, centered on the target particle 
position Xj. This grid covers the volume of the neigh- 
bor sphere. For each grid point, the normalized distance 
Xvoid to all neighboring particles can be calculated us- 
ing Equation ( 33 ) and the minimum value chosen. If the 
maximum value on the grid of x^oid > Cvoid, the position 
of the grid point with the maximum value is reported as 
a candidate void for particle creation. 

To speed up the calculation, we sieve the grid points 
lying within the neighbor sphere. We begin the search 
by initializing a large value on each grid point for the 
minimum value of Xyoid for that grid point. We then 
proceed by selecting each particle i in turn, and looping 
over all grid points. For each grid point, we calculate the 
normalized distance x^oid to the particle i. If its value 
for the grid point is less than the current minimum value 
on that point, we replace it with the newly calculated 
value for particle i. If the new value is less than c^oid, 
that grid point can be eliminated from the active list of 
candidates for void identification. We then move to the 
next particle and calculate its distance to the remaining 
active grid points, repeating the above procedure. After 
all particles have been sieved, if any grid points remain 
as void candidates, we report the one with the maximum 
Xvoid as a candidate for void creation. To hasten the 
operation, we first sieve the particles within A from the 
target position, then those within (3/2) A, and then the 
remaining particles, where A is the value for the target 
position. 

4.2. Clumps 

As the particles move, random fiuctuations will move 
them closer or farther from their neighbors. If two par- 
ticles approach each other too closely compared to A, 
they are essentially sampling the same field variable in- 
formation, and so are redundant. Because there are no 
restoring forces in the algorithm to separate nearby par- 
ticles, we instead remove any particle clumps of this sort, 
saving the computational cost of evolving the redundant 
particles. The question then remains of how to determine 
when a clump has formed. 

To do this, we define a scaled distance between two 
particles: 



{Xj Xj) 

XiXj 



(36) 



The nearest neighbor i to the target particle j is deter- 
mined. In turn its nearest neighbor is found. If they are 
mutual nearest neighbors and if 

they are candidates for deletion. We find that a value of 



^clurap 



= 0.12d„ 



(38) 



dition given by Equation ( 34 ). If it does we add a particle 



is suitable to prevent over-resolution. Among those two 
particles, we delete the one which was more recently cre- 
ated, retaining the particles with longer history to mini- 
mize the numerical diffusion from adaption. 

If particles are to be deleted, this is done so without 
considering whether that particle triggered the proposed 
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addition of a particle (that is, whether it is in a clump 
on the edge of a void). However, any proposed addi- 
tion resulting from the processing of that particle is still 
considered. 

4.3. Particle Creation and Deletion 

The first examination of all the active particles results 
in a proposed list of positions requiring particle addition, 
accompanied by the radius of the void detected. These 
proposals overlap, as each void may be detected by more 
than one particle. The list is exchanged by processes han- 
dling neighboring spatial domains, so that each process 
has a list of all the proposed additions within a distance 
of its boundary. The proposed addition list is then 
pruned, to select one position in which to add a particle 
within each void radius. To do this, each particle addi- 
tion proposal is compared to all other proposals within 
that spatial domain. If any other proposed location lies 
within its void radius, the values of the void radii are 
compared, and the proposal with the smaller void radius 
is rejected. 

We then create particles at the successfully proposed 
positions. Particles in this algorithm represent sample 
points, not discrete parcels of gas. When we add parti- 
cles, we are just sampling the continuous field variables 
at new positions. Therefore, considerations of conserva- 
tion do not enter this proc ess, unlike in particl e-splitting 
methods used in SPH (e.g. |Kitsionas fc Whitw orth 2002). 

The task of creating new particles needs to be load 
balanced among processors in order to handle situations 
where the memory required for new particles represents 
a large fraction of the total free memory in the parti- 
cle arrays. The new particles are then initialized in free 
spaces, on the processors to which they have been as- 
signed by the addition load balance procedure. As we 
have now deleted some particles, and added others to es- 
sentially random processors, a new load balance may be 
calculated among all particles, and the neighbor search 
data structure must be updated. Doing this on the entire 
particle list brings the new particles to optimal positions 
on the processors and provides neighbor information for 
the subsequent processing stages. 

New particles are initialized using a third order m ovin g 
least squares fit, as opposed to an interpolation (^3.1). 
This fit is centered on the position of the new particle. 

5. TIME STEPS 

The time step for each particle is set by taking the min- 
imum of five criteria. These are evaluated at the phase 
where new time derivatives are computed. The basic 
limit is the CFL condition for the stability of a forward- 
time-centered-space discretization. It is used here with- 
out explicit derivation as the general principle applies 
that the maximum stable time step must be short enough 
that a signal cannot cross a distance exceeding the local 
resolution A, 



At, 



CFL 



A 



CFL 



(39) 



where A^cfl is the CFL time step, Cqfl is the Courant 
number, which we usually take to be 0.3, Cg is the sound 
speed, and va is the Alfven speed. Note that, unlike 
an Eulerian code, the fiow velocity does not enter this 
equation. 



The use of the bulk viscosity fields Cs and Q require an 
appropriate time step constraint related to the diffusion 
term they introduce in the momentum equation. 



Ati; = 



Cc 



FL 



(40) 



Another time step constraint of the same form is applied 
based on the need for a sufficient number of time steps 
during a compression or expansion to allow for particle 
addition and deletion. This is the von Neumann time 
step 

where Cyc is a constant, which we usually take to be 2. 
The arbitrary form of the constant term Tr^Cyp comes 
from an analogy with the form of the von Neumann time 
step constraint by considering the von Neumann term 

C^vcA'(-V • V) + (-V • V) (42) 

where ()_)_ denotes that the expression is zero if the term 
contained is negative, as a diffusion operator and foUow- 
ing th e time step constraint from |Maron fc: Mac Low| 
(2009 Eq. 8). We also introduce a similar constraint 
based on the shear of the flow to allow for needed regu- 
larization, although the constraint on this vorticity time 
step much looser than in compression and expansion: 



At 



C, 



CFL 



VR 



107r2C2(.|V X 



v\ 



(43) 



The factor IOtt^Cvc is an ad-hoc scaling that in practice 
has been found to be sufficient. 
The time step limit assigned for a particle is 



At = min(AtcFL, At^, Atyc, AtvR)- 



(44) 



Each particle has an individually assigned time step. To 
tailor the algorithm to parallel implementation, a block 
time step scheme can be adopted. For such a time step 
scheme, the time steps actually used are rounded down 
to the nearest block time step interval. 

It is also necessary to ensure some degree of spatial co- 
herence to the time steps, so that disturbances propagate 
from short-time step particles to long-time step particles 
smoothly. At the end of the time update procedure for a 
particle, after the assignment of the new time step for a 
particle, if any of the neighbor particles has an end time 
greater than the target particle's new end time, a time 
step limit propagation procedure is triggered for the tar- 
get particle. The target particle's end time propagates 
to its neighbors, and if any neighbor's end time is fur- 
ther from the current time than twice the interval to the 
target particle's end time, then the neighbor's end time 
is set to this limit. 

6. ARTIFICIAL DIFFUSION 

Three types of artificial diffusion terms are used to sta- 
bilize the solutions to the equations modeled in Phurbas. 
These are terms acting as bulk viscosity, mass and inter- 
nal energy diffusion terms, and a term acting to diffuse 
magnetic monopoles. We find that these terms are suffi- 
cient to damp small scale fluctuations that would other- 
wise make the scheme unstable. 
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6.1. Bulk Viscosities 

We use two bulk viscosity fields. The first, Q, quenches 
small scale compressive motion in all areas of the flow, 
and evolves according to Equation The second, Cs, is 
a shock and discontinuity viscosity tn at ev olves according 
to Equation (10 1. Equations ^ and (10) each consist of 
three terms, a diffusion term, a source term, and a decay 
term. This configuration ensures that the bulk viscosity 
fields vary smoothly in time and space. 

The diffusion operator on the bulk viscosity fields 
= O.lSAcinax IS chosen to place a time step limit less 
stringent than the Courant limit from the hyperbolic part 
of the MHD equations. The bulk viscosity fields have 
source and decay terms and these terms have associated 
timescales. For the Q field the source and decay term 
timescales are the same, r/ = A/cmax- The Cs field is de- 
signed to fall off more slowly than it rises, which is partic- 
ularly advantageous calming post-shock oscillations, so 
rs_ = 20A/cinax and Ts+ = A/cmax- The use of a diffu- 
sion equation with source and decay terms to derive the 
artificial viscosity field here is analogous to the design in 
a disco ntinuous Galerkin method by [Barter fc Darmofal| 
( 2010 1 and the slow decay of the shock viscosity achieves 



a similar effect to the bulk v iscosity prescription used by 
Morris fc Monaghan[ ( |l997[ ). 
I'he Q source term is given by: 

S, =max(C^^(-a,F,)-HA', 
Ce\dM<^/p)\X'', 



Cpp-'\d.Aip)\>^'^/c[+^a, 
CpP''\dMP)\^'Vc[+^). (45) 

The first term has the form of the conventional von Neu- 
mann artificial viscosity, with Cvn = 2. The second 
term responds to changes in the specific internal e nergy, 
in a rn anner similar to the dissipation introduced by |Price 
(2008). This form is a trigger on the size of the second 
derivative of the specific internal energy, and Ce = 0.1. 
The third term is constructed analogously to the second, 
but using the Laplacian of density and Cp = 1.0. The fi- 
nal term is again constructed analogously to the second, 
but using the Laplacian of pressure and Cp = 3.0. The 
constants Cvn, Ce, and Cp can be tuned for a partic- 
ular problem, with smaller values being preferable, but 
the values given here have proven to be sufficient for most 
problems. 

The mass and internal energy diffusion terms are cou- 
pled to the Cs field, with a strength set by the constants 
Hp = Ha = 5 X 10~^. For stability, it would be preferable 
to have a small scale mass and internal energy diffusion 
(such as a hyperdiffusion) active everywhere in the flow, 
but we have not found a formulation of such a term that 
is sufficiently accurate to yield reasonable mass conser- 
vation results. 

6.2. Magnetic Divergence Diffusion 

As Phurbas solves the equations of MHD, the issue of 
magnetic monopole errors must be treated. The primary 
problem caused by monopole errors in schemes of this 
type is numerical instability. The interpolating moving 
least squares derivative estimates may return derivatives 
of the magnetic field that do not satisfy V • B = 0. Over 



several time update cycles, these estimates may lead to 
the local creation of a net magnetic monopole character 
to the field. In practice we have found that a diffusive (or 
parabolic) correction is sufficient to prevent the growth 
of this monopole character of the magnetic field (see tests 
in Paper II). 

For each particle, a V-B field is defined. The value is 
simply reset each time step to the value oiV -B derived 
from the fit to the magnetic field. The derivatives of 
the V • B field derived from the fits are then used to dif- 
fuse V • B, generally resulting in a reduction of its value. 
These fitted values and derivatives of V • B are less noisy 
than values and derivatives of V • B derived directly from 
fits to the magnetic field. 

The diffusion term for particle j is 

V(V • B,), (46) 

where r]max is the maximum diffusion coefficient possible 
under the stability criterion 



At < 



(47) 



from Maron fc Mac Low 
by Equation 



2009[ Eq. 8). The term given 
45[) is added to the right hand side of the 



induction equation (Equation [2]). in the first time deriva- 
tives used in the second-order predictor-corrector scheme 
for the evolution of the magnetic field. The effect of this 
is that the V-B diffusion operator is integrated with a 
first-order predictor-corrector scheme. The V B diffu- 
sion rjjnax is computed each time the fields are fit, which 
occurs at times that are the end of one time step and the 
beginning of the next. The time step used to define rjmax 
is the time step that has its end at the instant r]rnax is 
calculated, i.e. the previous time step. Thus, the rjmax 
used in the predictor stage of the time integration of a 
particular step is different from the rj^^ax later used in 
the corrector stage of the same time step. This diffu- 
sion is not conservative, but the Phurbas discretization 
only preserves the conservation in the MHD equations to 
truncation error levels, and the V-B operated on by this 
diffusion is, by definition, only created below truncation 
error levels. We find that since the canonical form of the 
Lagrangian MHD equations that we use treats V • B as a 
passively advected scalar, the presence of small amounts 
of V • B does not destabilize the solution. 

7. SUMMARY AND DISCUSSION 

7.1. Summary of the Algorithm 

We now summarize the conceptual steps of the algo- 
rithm. The operations described here are actually often 
broken into multiple phases to enable efficient paralleliza- 
tion. Future implementers of the algorithm should con- 
sider the specific needs of each operation when designing 
data structures and communication patterns. 

• Build the neighbor-finding data structure, such as 
a particle tree. 

• Balance particles among processors. 

• Identify target particles for evolution. 

• Use the tree to assemble all neighbors within a ra- 
dius r/ of the target particles. 
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• Check for voids. Com plete particle addition for 
qualifying voids (§ 4.1 1. 



Check if there are mutual near est n eighbor pairs 
that are too close. Delete one (§ |4.2[). 



Evaluate the polynomial fit weights (§ 3.2) 



• Compute the local polynomial fit to derive values 
and sp atial derivatives at the location of each par- 
ticle (§[3l|). 



• Use the polynomial coefficients to evaluate the 
MHD equations for the Lagrangian time derivatives 
including diffusivity terms (§l2j). 

• Use the polynomial coefficients to evaluate the gov- 
erning equations of the bulk viscosity fields Q and 
Q for the Lagrangian time derivatives (§[2]). 

• Use the ti me d erivatives to correct the previous 
time step (§ 3.3). 



Evaluate the resolution scale A for each particle 
position (§[2]). 



Evaluate the size of the next time step (§[5]). 

Use the time derivatives to p redict forward in 
to the next time step (§|3.3). 

Restrict local time step variations (§[5]). 



7.2. Effective Resolution 

To understand the effect of varying the effective reso- 
lution parameter A on the numerical resolution, consider 
a one-dimensional uniform grid with a grid spacing of 
unity, initialized with a field variable having values given 
by the Fourier mode sin(7rfca;). The maximum wave num- 
ber k of this mode that can be expressed on this grid is 
fc = 1, the Nyquist wave number. In order to be able 
to calculate realistic derivatives with finite differences, k 
must be le ss than unity , and th e pre cision increases as k 



decrea ses. Maron et al. (2008) and Maron & Mac Low 



( 2009 ) evaluate the effective precision of hnite ditterence 



schemes with varying stencil sizes. They find that for a 
stencil radius of {1,2,3,4}, finite differences can be cal- 
culated with a relative precision of ~ 1 percent up to 
a wave number of A; ~ {1/8,1/4,2/5,1/2}. Given that 
derivatives are more easily calculated on a grid than for 
irregular particles, we take this as an upper limit for what 
we can expect from particles. Since a 3D, third-order, 
polynomial corresponds to a 5-point (or stencil radius 
2) ID finite-difference scheme, we expect fc ~ 1/4 to be 
the limit of resolution, corresponding to a wavelength of 
- 8A. 

8. CONCLUSION 

We have described Phurbas, an adaptive, Lagrangian, 
meshless algorithm for MHD. The algorithm is described 
for the specific case of the MHD equations, but can be 
easily generalized to other hyperbolic systems, as the fit- 
ting, time integration, and stabilization procedures do 
not rely on particular properties of the MHD equations. 
The central principle of the algorithm is that the solu- 
tion and its spatial derivatives are derived from a high- 
order, interpolating, polynomial fit to a set of particles 



that are merely Lagrangian sample points in the flow, 
not mass elements as in SPH or finite volume methods. 
This allows for significant flexibility in the design of the 
algorithm, and the implementation of additional physical 
processes. Particle addition and deletion is required to 
prevent the growth of voids or clumps. This naturally al- 
lows the numerical resolution to be fully adaptive based 
on user specified criteria. The Lagrangian nature of the 
code means that particles can be evolved with time steps 
dependent only on the nature of the local flow, and that 
numerical diffusion is Galilean invariant. The version 
described here is just one subset of the many available 
options. Paper II describes a parallel implementation 
and tests of Phurbas that demonstrate accuracy compa- 
rable to that of third order grid codes on subsonic and 
supersonic problems. 

A few theoretically desirable improvements to the 
scheme can already be identified. Though the Q field 
is sufficient to modify the equations to give reasonable 
results in tests, it would be preferable to only use stabi- 
lizing viscosities that scale as A^ or a higher power (hy- 
perviscosities) to give faster convergence of the modeled 
equations to the limit of ideal MHD. Possible avenues 
though which such an effect could be achieved are use of a 
different interpolation scheme than moving least squares, 
and/or a time or spatially dependent version of Q so that 
it couples only to the shortest wa velength or shortest 
timescale motions. Kuhnert ( 1999 ) showed methods for 
introducing upwindmg into a scheme similar to Phurbas, 
a modification that may reduce the need for Q. Addi- 
tionally, in defining we have used a simple formulation 
where the effects decay on the same timescale and re- 
flect the effects of discontinuities (shoc ks and otherw ise) 
through the same parameter. In SPH, iPricel p012 ) and 
,Rosswog, (2009| ) have found that separate parameters for 
each discontinuity are useful. Analogously in the frame- 
work here, (s could be broken into three flelds, though 
this would come at some cost. 

The second order Hermite predictor-corrector scheme 
is particularly useful as it only requires computations of 
spatial derivatives at the beginning of each time step, and 
the predictor half of the integration can be done without 
knowing the end-time of the time step. It however has 
the drawback that first and second time derivatives of 
the field variables must be obtained. In general, these 
analytic expressions could be very complicated. A time 
integration scheme, such as a Runge-Kutta method, that 
uses only first time derivatives may be preferable in this 
sense. 
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APPENDIX 
STABILITY OF A MODEL SCHEME 



To determine the stability requirements for our meshless method, we here present a stability analysis for a simple 
model that captures the essential points of our algorithm. We analyze a numerical scheme for approximating the 
solution of the diffusion equation 



du 
~dt 



D 



(Al) 



We discretize u{x) on a grid with spacing Ax. To find derivatives of u we use a moving least squares approximation. 
Shifting the origin to grid point xq the polynomial approximation of u{x) is given by 



U{x) = apX^ . 

p=0 



(A2) 



The coefficients Op are determined by minimizing the sum of the square of errors E2 at 2A^ + 1 neighboring grid points, 
given by 



N 
j = -N 

Combining these two expressions, and using the definition of the grid point positions xj — jAx gives 

N P 
j = -N p=0 

We can find the the minimum of E2 by setting dE2 jdaq — yielding 

N N P 

J2 fAx'^u,^ J2 Y^P^"^'^'''^''' 



(A3) 



(A4) 



(A5) 



j=-N 



j = -Np=0 



which is a system of equatio ns fo r q = 0..P that can be solved for the polynomial coefficients Up. The second derivative 
of the polynomial Equation (A2) at a; = is 2a2. So, we can write a scheme to update the solution to Equation ( Al I 
as 



,n+l 



flo + 2Da2At, 



(A6) 



where uq = u{xo)- In each step we replace the value Mq with the fit qq and use the second deriv ative of the moving least 
squares fit to construct a forward-Euler type time update. This is a Maron & Howes ( |2003 ) type scheme for solving 
the diffusion equation. As an example, for third order polynomials (p = 3) and iV — 4 this scheme is specifically 



,n+l 



Uo 



f 59 10 DAt\ 59 A , , 5 A 2 . 



DAt 



V231 231 A.t2 ) 

4 



231 



-1 
231 



(A7) 
(A8) 



To perform a von Neumann stability analysis, we substitute = ^"e 



71 Ak^\x 



and solve for ^ 




59-10 



,DAt 
A^ 



^(59- 5j2)cos(jfcA2;) 
i=i 



DAt 
A^ 



E 



-1 f_ 

231 154 



cos(jfcAx). 



(A9) 



We can evaluate this expression numerically, plotting the amplification factor |^| for each wave number kAx and time 

(2003) type schemes using a 



Maron & Howes 



step This is shown in Figure The fundamental trouble with 

moving least squares fit is the region at low wavenumbers and small time steps where the rnagnitiide of the amplification 
factor is greater than unity. In this region the scheme is unstable. 
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Figure 1. Von Neumann stability analysis of model schemes, showing the amplification factor \^\ as a function of perturbation wavelength 
k and time step At. bo th appropriately normalized to the grid. The unstable region with values of > 1 is shown in white. Left: the 
[Maron fc Howes] | |2003| type scheme for the diffusion equation. Note the instability of low wavenumber perturbations kAx < 0.4. Middle: 
Lieast squares mterpolant scheme for the diffusion equation. Right: Least squares interpolant scheme for the advection-diffusion equation. 
Both the latter schemes show a region of stability at small enough time step for all wavenumber perturbations. 



If instead we use a moving least squares interpolant instead of just a fit, the system of equations for the coefficients 
is given by 



JV N P 

3 = -N j=-Np=l 



(AlO) 



for q = 1,2,3. Here we have eliminated the coefficient ao by forcing the moving least squares approximation to 
interpolate uq. Then we can write a forward-Euler type scheme using the second derivative of this interpolant as 



71+1 n 



+ 2Da2At. 



For p — 3 and = 4 this gives the scheme 



,n+l 



Uq ■ 



DAt v;^ 

3=0 



2uo) 



Uq 1 



eODAt 



DAt 



Again, performing a von Neumann stability analysis, we substitute u" — ^"e and solve for ^ 

4 



1 



eODAt 



DAt 



354Aa;2 / 354Aa;2 



E j^2cos(jfcAx). 



(All) 

(A12) 
(A13) 

(A14) 



i=i 



Th e magnitude of £ in this case is shown in the center panel of Figure (fTl) for the case D = cAx. In stark contrast to 
the Maron & Howes (2003) type scheme constructed with a non-interpolating, moving least squares fit, a significant 
region of stability (K| < 1) exists at small timestep for all wavenumbers kAx. 
For the advection equation 



du du ^ d u 
dt dx dx-^ 



we can write the scheme 



n+l n 



+ cAtai + 2ca2AxAt. 



(A15) 



(A16) 



Here we have set the diffusion parameter to be scaled by the grid resolution Ax. With the interpolating moving least 



14 



squares coefficients ai and 02 as in the diffusion problem above we have 



4 



4 



A von Neumann stabihty analysis yields 



cAxAt -r-^ .? / N / » 1 „N 



4 4 



4 



^j22cosO-fcAx). (A18) 



354Aa; 



Note that the second term is imaginary, arising from the advection operator, so if the contributions from the diffusion 
operator, the part of the first term, and the final term, were dropped, the scheme for the pure advection problem 
would be unconditionally unstable. We plot the magnitude of the amplification factor |^| in Figure ([l]). The addition 
of the diffusion operator has stabilized the advection problem for sufficiently small time steps cAt/Ax. 



TIME INTEGRATION 



Time integration proceeds using a system of Hermite predictor-cor rector formulas for field va riable and position 
updates. This scheme is a lower order version of that presented in Nitadori & Makino| (|2008|). As it is has not 



previously been described in the literature, we present here a brief derivation of the integration method. 

We predict field values Qp^i+i at time t+At using a Taylor series expansion around their values at time t incorporating 
their time derivatives calculated after the prediction phase of the previous time step, giving 

Qp,i+i = qc,i + DtqpsAt + ]^DttqpAAt^ (Bl) 

where (7c,i is the corrected field value from the previous time step at time t. Af ter calculating the values of g^.i+i, 
we use them to calculate Dtqp^i+i and Dttqp,i+i, as given by Equations ([5)-(10l, and those in Appendix [C| These 



derivatives will be used in the corrector stage of the current time step and in the predictor stage of the next time 
step. The stabilizing diffusion terms, linked to the (s and Q fields are proportional to the resolution parameter A. As 
asymptotically, the time step varies as A itself due to the CFL limit, we only integrate these terms to first order in 
time, so the combined space-time error is second order in A. 

By using higher time derivatives, a Hermite scheme of time integration depending only on field variable values from 
the previous time step can be constructed. This saves storage, avoids complex start up procedures, and simplifies the 
use of individual particle time steps in comparison to predictor-corrector schemes based on Newton interpolation (such 
as Aarseth and Adams-Bashforth-Moulton schemes) that require storage of field values from earlier time steps. The 
Hermite corrector stage is constructed as 

qc,i+i = qc.i + / fc{T)dT. (B2) 
Jo 

We choose the function /c(t) to be a Hermite interpolation, that is, a polynomial that interpolates Dtq and Dttq 
at each end of the time step. We further choose to simplify the formalism by designing the polynomial to be time 
symmetric about t + At/2, so that 

AA „ / AA^ , / AA^ 



fcir) = fo + h[r-—j+.f2[T-—j +f,l^r-—j . (B3) 

We determine the coefficients /o, /i, /2, /a by using four constraints: at r = 0, fdT) must have a value of Dtqpj, and 
a time derivative Dttqp,i', while at r = Ai the valu e and the derivative must be Dtqp^i+i and Dttqp,i+i, respectively. 



However, evaluating the integral in Equation (B2), the time symmetry we chose yields the simple result that th e /i 



and /3 terms integrate to zero regardless of the values of their coefficients. Performing the integral in Equation (B2| 
yields a correction stage 

<7c,i+i = qc,i + ^{Dtqp,i + Dtqp^i+i)At + ^{D^qp^i - At9p,i+i)Ai^. (B4) 

Velocity V is treated as an independent set of field variables, so for particle positions x there are three time derivatives 
of information available, as well as corrected values of velocity from the beginning and end of the current time step. 
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Therefore, for the predictor stage for the position we use a third order Tayfor series incorporating the best information 
available at time t, 

Xp^^+i = x^,, + Vc,^At + ^AT^p,^A^^ + ^DttVp^^M^. (B5) 

Similarly to the field variable integration we choose a time-symmetric, Hermite interpolating function, so that the 
corrector stage is 

Xc.i+l = Xc,'j. + / gc{T)dT. (B6) 

The function gc{T) is again a polynomial centered on i + At/2, that now interpolates through Vc.i, DtVp^i, and DttVp,i 
at T = 0, and through Vc,i+i, DtVp^i+i, and DttVp^i+i at r = At. (The availability of Vca+i occurs because we have, 
at this point, already updated the field variables.) Applying these constraints allows us to evaluate the coefficients in 
the interpolating polynomial 

9c{r)^go+9i{T~ —)+g2{T- —f +g^{T- —f +gi{T- —f + g^{T~ —f. (B7) 



Evaluating the integral in Equation (B6), the gi, g^, and 55 terms are zero due to the choice of time symmetry, and 
so the correction stage for particle positions is 

Xc,^+l = x,^, + + V,,,+i)At+ ^{DtVp,, ~ DtVp^,+i)At^ + ^{DuVp^, + DuVp,^+i)At^ . (B8) 

It can be useful to split the integration in to a background fiow and perturbations, for example in computing the 
dynamics of a steady-state cylindrical flow. In this case we define the perturbation velocity field as V' = V—^r4> where 
r is the two-dimensional radius from the center of the cylinder, and ^{r) is the angular velocity of the background 
flow. We denote the components of the circular radius to the point {xc.i,yc,i) as rc,i,x-,'fc,i,y, and the angular velocity 
of the background flow at this radius is Qp. Then, the predictor step with the background flow separated from the 
perturbation velocity V' is 

Xp,;+i = r,^,,x cosiflpAt) - r,,,,y sin(r!pAt) + V^.^^At + ^DtV^^^ ^^At^ + ^DuV;^,,,At^ (B9) 

yp,^+l = rc,^,y cos(r!pA<) + r,,,,, sin(r!pAt) + K'z,„Ai -f ^^DtV^^^^yA^ + ^DuV^^.^yAt' . (BIO) 

We then add the background flow state back into the field variables used in Phurbas to calculate time derivatives for all 
fields in the usual inertial reference frame with y not Y . Then, to transform the time derivatives to time derivatives 
of the perturbation velocity we use 





-l.x 


= DtVp^i+i^x ~ 


^p'^p,i+\^x 


(Bll) 




n,y 


= DtVp,i^i^y - 




(B12) 






~ DttVp^iJ^i,x 


^p^P,^+l,y 


(B13) 




-l,x 


— DttVp^i-\-i,x 


^p^p^i+l.x- 


(B14) 



The perturbation velocity V can be integrated directly using these perturbation time derivatives. Using the normal 
corrector for the field variables and the perturbation velocity, the corrector step for position is then 

Xc.+i = rc,.,. cos(f]pAt) - r,,,^y sin{%At) + \{Vi^,^x + K«+i.JA^ + ^(AX^p',.,, - AV^'.^+i.JAt^ 

+ Y^(AtV;',,,. + AtV;;,+i,JAi3 (B15) 

yc,z+i = rc,..j, cos(r!pAi) -t- r,,,,, sin(17pAt) + \{yl,^y + VU^x,y)At + ^{DtVl^^.^y - Al^p',,+i,y)Ai2 

+ -f DuVi,^r^y)At\ (B16) 

To prepare for the next step, it is necessary to shift the perturbation velocity and time derivatives into the correct frame 
that will be used for the next step. The new frame at the corrected position (xc,i+i, yc,i+i) has radius components 
(rc,i+i,a;, rc,i+i,y) and the angular velocity of the background fiow at this radius is fie- The transformations made to 
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these quantities are: 







(B17) 






(B18) 






(B19) 






(B20) 






(B21) 






(B22) 



This shifts V' into the accelerating reference frame used for the following predictor step. 

SECOND TIME DERIVATIVES OF MHD EQUATIONS 



For the second order predictor-corrector scheme § 3.3 we need both the first and second Lagrangian time derivatives 
of the MHD equations Q-Q. This appendix gives formulas for the required second time derivatives. 
The form for a Lagrangian time derivative Dt in terms of partial derivatives d is; 

Dtdq = dtdq + • Vg (CI) 
= dtdq + • Vg] - dVAq (C2) 

= dtd.q + d, ~~ {d,v,){d,q) (C3) 

= dDtq-[dV,){d,q) (C4) 

Applying this to the MHD equations ([l])-(|4]) gives the second time derivatives needed for the second order predictor 
corrector scheme. We start with the velocity equation (Equation [ij . Taking its Lagrangian time derivative, the first 
term on the right hand side becomes 

Dti^-p-^d,P) = -p-\d,DtP - djVadaP) + p-^Dtp)djP. (C5) 

Inserting this, the second derivative of velocity is 

DttVj = [djDtP - d.VAP) + p-^d,P)DtP + Dt{p-\e,abeacdidcBd)Bb)). (C6) 

This equation can be further reduced. The two pressure-dependent terms depend on the equation of state. For a 
gamma-law equation of state, P = (7 — l)a, 

DtP = (7 - l)Dt(7 = (7 - l)(-(a + P)d,Vi), (C7) 

and so 

d.DtP = (7 - 1) (-(9,^ + d,P)daVa + P)dajVa) . (C8) 

For an isothermal equation of state P ~ c^p, 

DtP^clDtp^cK-pd^Vi), (C9) 

and so 

d.DtP = cl {-d.pdaVa - pdajVa) . (CIO) 

The magnetic term reduces to 

Dt{p^^{ejab£acd{dcBd)Bt)) ^ {Sjab£acd{dcBd)Bb){~p~'^Dtp) 

+ {SjabSacd imB,)deVd + B,d,,Vd - {d,Bd)d,Ve - SdSceK - {d,V,)d,Bd) Bb 

+{d,Bd)DtBb)) , (Cll) 
while the Lagrangian time derivative of the gravitational force 

Dt{+Gj) = dtG, + VAGj. (C12) 
Taking the Lagrangian time derivative of the induction equation (Equation [2]) gives 

DttBj ={DtBi)d,V, + B,{d,DtV, - diVAVj) 

- (DtBj)d^V^ - Bj{d,DtV, - d.VkduVi). (C13) 
The Lagrangian time derivative of the internal energy equation is 

Dtta = -{a + P){d,DtV, - d^V.djV,) - (Act + DtP)d,V,, (C14) 
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where the required pressure-dependent expressions have aheady appeared above. If we have a barotropic or isothermal 
equation of state then this equation is not used, of course. Finally, the Lagrangian derivative of the continuity equation 
(Equation |4| is 

DttP = -{Dtp){d^Vi) - p{d,DtV, - d.VjdjV,). (C15) 
The second term on the right hand side can be expanded as 

^^DtV, =e,abeacd {-p-\d,Ba)Bi,d,p + p-^dMB^, + d.Bad.B^) 

+ p-^d,Pd,p-p'^duP + d,G,. (C16) 



